##########################################################################################

library(data.table)
library(optparse)
library(parallel)
library(dplyr)
library(RColorBrewer)
library(immunedeconv)

##########################################################################################

option_list <- list(
    make_option(c("--cibersort_path"), type = "character"),
    make_option(c("--rsem_file"), type = "character"),
    make_option(c("--gtf_file"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    cibersort_path <- "~/tools/StandTools/Cibersort/"
    rsem_file <- "~/20220915_gastric_multiple/rna_batch1/analysis/RSEM/CombineTPM.tsv"
    gtf_file <- "~/ref/GTF/gencode.v19.ensg_genename.txt"
    out_path <- "~/20220915_gastric_multiple/rna_batch1/analysis/images/cibersort"

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

cibersort_path <- opt$cibersort_path
rsem_file <- opt$rsem_file
gtf_file <- opt$gtf_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)

##########################################################################################
## https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5895181/
## 文章提出Standard RNA-Seq expression quantification metrics, 
## such as fragments per kilobase per million (FPKM) and transcripts per kilobase million (TPM), 
## are suitable for use with CIBERSORT.



##########################################################################################

dat_tpm <- data.frame(fread(rsem_file))
colnames(dat_tpm) <- gsub( "X" , "" , colnames(dat_tpm))
dat_gtf <- data.frame(fread(gtf_file , header = F))
colnames(dat_gtf) <- c("gene_id" , "Hugo_Symbol")

dat_tpm$gene_id <- sapply( strsplit(dat_tpm$gene_id , "[.]") , "[" , 1 )
dat_tpm_all_out <- merge( dat_gtf , dat_tpm)

#去除重名基因
dat_tpm_all_out <- dat_tpm_all_out[!duplicated(dat_tpm_all_out$Hugo_Symbol),]
dat_tpm_all_out <- dat_tpm_all_out[,-1]
colnames(dat_tpm_all_out)[1] <- "gene_id"

##########################################################################################

cibersort_r <- paste0(cibersort_path , "/Cibersort.R")
lm22 <- paste0(cibersort_path , "/LM22.use.txt")

source(cibersort_r)
sig_matrix <- lm22

## 只提取免疫浸润用到的基因
dat_lm22 <- data.frame(fread(lm22))
dat_tpm_all_out <- subset(dat_tpm_all_out , gene_id %in% dat_lm22$Gene.symbol)

mixture_file <- paste0(out_path , "/CombineTPM.Hugo_Symbol.tsv")
write.table(dat_tpm_all_out , mixture_file , row.names = F , sep = "\t" , quote = F)


##########################################################################################
## https://liulab-dfci.github.io/RIMA/Infiltration.html#cibersort

#Run CIBERSORT abs 
#The number of permutation
cibersort_perm = 1000
#Quantile normalization of input mixture, default = FALSE for RNA-Seq data
cibersort_qn = FALSE
#whether to apply absolute mode in cibersort
cibersort_abs = TRUE
#sig.score = for each mixture sample, define S as the median expression,level of all genes in the signature matrix divided by the median expression level of all genes in the mixture. Multiple cell subset fractions by S.
cibersort_abs_method = "sig.score"
res_ciber <- CIBERSORT(sig_matrix = sig_matrix, mixture_file = mixture_file, perm = cibersort_perm, QN = cibersort_qn)

##########################################################################################

#res_ciber <- CIBERSORT(sig_matrix, input, perm = cibersort_perm, QN = cibersort_qn, absolute = cibersort_abs,
#                       abs_method = cibersort_abs_method)

out_name <- paste0(out_path , "/cibersort.tsv")
result <- data.frame(res_ciber)
result$Sample <- row.names(res_ciber)
result <- result[,c(ncol(result),1:ncol(result)-1)]
write.table( result , out_name , row.names = F , sep = "\t" , quote = F )
